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Figures 


1.  Discontinuities  in  the  mantle  transition  zone  at  410-  and  660-km  depths  cause  a  pair  of 
triplications  in  P-  and  S-wavc  arrivals  between  about  15  and  30  degrees.  The  ray  paths 
(top)  are  color-coded,  and  correspond  to  the  different  branches  in  the  travel-time  curve 
plotted  below  (reduced  at  10  s/degree).  The  AB  branch  consists  of  direct  waves  that 
bottom  above  the  410-km  discontinuity  (red  solid).  The  BC  branch  reflects  at  410  km 
(red  dashed).  The  CD  and  DE  branches  bottom  above  and  reflect  off  the  660-km 
discontinuity  (blue  solid  and  dashed).  The  EF  branch  bottoms  in  the  lower  mantle.  This 
figure  is  adopted  from  Shearer  (2000).  1 


2.  Lateral  and  vertical  variation  in  seismic  velocity  structure  (Masters  et  ah,  2000)  and 

mantle  topography  (Flanagan  and  Shearer,  1998)  cause  significant  variation  in  the  upper 
mantle  triplications.  The  location  of  this  profile  is  shown  in  the  map  at  the  center. 
Topography  is  exaggerated  by  3  times.  2 

3.  Distance-dependent  Earth-response  terms  for  R-wave  amplitude  from  10  to  90 

degrees.  4 

4.  A  graphical  representation  of  the  sources  of  seismic  signals  (event,  Earth  response 

tenn,  receiver  tenn,  and  noise)  that  we  attempt  to  decompose.  5 

5.  Examples  of  three  different  types  of  perturbations  to  transition  zone  structure:  (a) 

transition  zone  thickness,  (b)  velocity  contrast,  and  (c)  interface  thickness.  The 
corresponding  WKBJ  maximum  amplitudes  for  P  waves  (d-f)  change  as  a  function  of 
distance.  The  models  are  shown  as  black,  blue,  and  red  lines.  The  data  are  the  dashed 
lines.  5 

6.  Parameterization  of  the  P  velocity  model.  Each  interface  (410,  520,  &  660)  can  vary  in 

depth,  thickness,  and  velocity  contrast.  7 

7.  The  (red)  best  fit  model  after  100  generations  of  an  NGA  simulated  inversion  yields  a 

similar  result  to  a  similar  parameter  search  based  on  (blue)  PdP,  SdS,  Pds,  and  Ppdp 
reflected  phases  (Lawrence  and  Shearer,  2006).  7 

8.  The  time-  and  distance-dependent  Earth-response  function  calculated  from  all  FARM 

P- wave  data  using  the  spectral  deconvolution  method.  The  box  outlines  the  region  of 
triplicated  data.  Blue  is  positive  amplitude,  red  is  negative.  Various  phases  are  indicated 
as  P,  PP,  and  PcP.  9 

9.  The  time-  and  distance-dependent  Earth-response  function  calculated  from  all  FARM 

R-wave  data  using  the  iterative  deconvolution  method.  The  box  outlines  the  region  of 
triplicated  data.  Blue  is  positive  amplitude,  red  is  negative.  Various  phases  are  indicated 
as  P,  PP,  and  PcP.  1 0 
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10.  The  amplitudes  of  the  globally  stacked  Earth  response  function  using  iterative 

deconvolution  at  time  zero  for  each  distance  (yellow)  grossly  match  the  average  estimates 
from  simple  maximum  amplitudes  for  E-wavcs  (red).  The  black  dots  illustrate  the  scatter 
in  the  P- wave  amplitudes.  1 1 

1 1 .  This  figure  shows  (a)  the  synthetic  Earth-response  function  for  the  seismic  velocity 

model  found  by  Lawrence  and  Shearer  (2006)  (LS06),  and  deviations  in  Earth-response 
functions  relative  to  LS06  as  a  result  of  (b)  a  deeper  660  interface,  (c)  a  larger  Vp^o 
contrast,  and  (d)  a  thick  410  interface.  Blue  is  positive,  red  is  negative.  12 

12.  A  comparison  showing  that  the  Earth-response  function  calculated  from  the  model  of 

Lawrence  and  Shearer,  (2006)  is  very  similar  to  the  global  Earth  response  function.  These 
figures  were  calculated  with  decomposition  using  spectral  division.  13 

13.  The  Earth  response  tenn  for  a  single  event,  as  (a)  traces,  and  (b)  as  the  zero  time 

amplitude  (x’s)  as  a  function  of  event-to-station  distance.  14 

14.  The  relative  scaling  factor  between  the  zero-time  Earth-response  amplitudes  of  a 

single  event  and  the  zero-time  Earth-response  amplitudes  of  the  global  dataset  plotted 
versus  body  wave  magnitude.  1 5 
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1.  SUMMARY 


This  project  studies  the  effects  of  heterogeneous  transition-zone  structure  on  seismic 
travel  times  and  amplitudes.  Accurate  predictions  of  P-  and  S- wave  travel  times  and 
amplitudes  at  distances  between  13  and  30  degrees  are  hindered  by  the  sensitivity  of  the 
multiple  travel-time  branches  at  these  distances  to  variable  structure  in  the  mantle 
transition  zone.  Both  discontinuity  topography  and  bulk  seismic  velocity  anomalies 
perturb  seismic  ray  paths,  which  causes  focusing  and  defocusing  effects  on  wave 
amplitudes  as  well  as  travel-time  anomalies.  However,  travel-time  and  amplitude 
information  is  critical  for  locating  and  estimating  magnitudes  of  target  events.  By 
comparing  regional  variations  of  triplication  amplitudes  and  travel  times  with  predictions 
of  3D  seismic  velocity  models,  it  is  possible  to  obtain  improvements  in  mantle  transition- 
zone  models,  as  well  as  in  the  estimated  locations  and  magnitudes  of  recorded  events. 
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Figure  1.  Discontinuities  in  the  mantle  transition  zone  at  410-  and  660-km  depths  cause  a 
pair  of  triplications  in  P-  and  5-wave  arrivals  between  about  15  and  30  degrees.  The  ray 
paths  (top)  are  color-coded,  and  correspond  to  the  different  branches  in  the  travel-time 
curve  plotted  below  (reduced  at  10  s/degree).  The  AB  branch  consists  of  direct  waves 
that  bottom  above  the  410-km  discontinuity  (red  solid).  The  BC  branch  reflects  at  410 
km  (red  dashed).  The  CD  and  DE  branches  bottom  above  and  reflect  off  the  660-km 
discontinuity  (blue  solid  and  dashed).  The  EF  branch  bottoms  in  the  lower  mantle.  This 
figure  is  adopted  from  Shearer  (2000). 
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The  effect  of  the  mantle  discontinuities  at  410-  and  660-km  depth  is  shown  in  Figure 
1,  which  depicts  a  series  of  P- wave  ray  paths  between  13  and  35  degrees.  Retrograde 
branches  result  from  the  velocity  jump  at  each  of  the  discontinuities,  causing  the  familiar 
double  triplication  centered  at  about  20  degrees.  However,  the  positions  of  these 
branches  are  very  sensitive  to  the  discontinuity  depths,  which  typically  vary  by  30  km  or 
more  (e.g.,  Shearer,  1990,  1991,  1993;  Shearer  and  Masters,  1992;  Flanagan  and  Shearer, 
1998a,b;  Gu  et  al.,  1998),  and  to  the  size  of  the  velocity  jumps  at  the  discontinuities, 
which  also  exhibit  considerable  variation  (e.g.,  Melbourne  and  Helmberger,  1998; 
Shearer  and  Flanagan,  1999;  Chambers  et  al.,  2005).  In  addition,  there  is  some 
sensitivity  to  3D  velocity  variations  (e.g.,  Figure  2)  in  the  mantle  transition  zone. 
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Figure  2:  Lateral  and  vertical  variation  in  seismic  velocity  structure  (Masters  et  al.,  2000) 
and  mantle  topography  (Flanagan  and  Shearer,  1998a,b)  cause  significant  variation  in  the 
upper  mantle  triplications.  The  location  of  this  profile  is  shown  in  the  map  at  the  center. 
Topography  is  exaggerated  by  3  times. 

We  compare  observed  P-wave  travel  times  and  amplitudes  with  predictions  based  on 
the  best  current  models.  These  tests  illuminate  the  strengths  and  weaknesses  of  these 
models,  from  which  we  determine  if  the  models  are  sufficient  to  explain  anomalous  travel 
times  and  amplitudes  at  these  distances  (and  consequently  check  to  see  how  robust  our 
current  event  locations  and  amplitude  predictions  are).  This  project  compares  times  and 
amplitudes  recorded  at  both  temporary  and  permanent  seismic  stations  for  all  triplication 
phases,  not  just  the  first  arrivals.  The  triplications  cause  difficulties  for  traditional 
methods  of  source  location,  magnitude  estimation,  and  inversions  for  3D  velocity 
structure.  Therefore,  triplication  data  are  typically  avoided  for  seismic  tomography  and 
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source  location.  However,  details  of  the  triplications  can  be  useful  for  characterizing 
heterogeneity  and  calculating  better  source  locations.  With  sufficient  data  the  secondary 
branches  can  be  used  to  model  variations  in  discontinuity  topography  and  seismic 
velocity.  Once  the  structural  properties  are  detennined,  the  high  sensitivity  of  triplication 
travel-time  and  amplitude  fluctuations  to  source  location  and  magnitude  make 
triplications  ideal  for  determining  source  characteristics.  However,  without  an 
appropriate  transition  zone  model,  these  times  cannot  yield  accurate  times  or  amplitudes. 

The  relevance  of  our  results  for  nuclear  test  monitoring  is  that  better  models  of 
transition-zone  structure  will  reduce  source  location  and  magnitude  uncertainties.  The 
anomalous  arrival  times  and  amplitudes  between  13  and  30  degrees  currently  limit  the 
usefulness  of  regional  phase  data  in  calculating  accurate  source  locations  and  event 
magnitudes.  However,  records  from  closer  distances  are  not  likely  to  be  available  in 
many  parts  of  the  world  and  small  magnitude  events  are  often  not  well  recorded  at  longer 
distances  because  of  the  sharp  drop  in  P  and  S  amplitudes  that  occurs  just  beyond  30 
degrees.  Thus  unraveling  the  complexities  of  the  travel-time  triplications  and  improving 
our  models  of  the  transition  zone  are  likely  to  be  critical  for  accurate  monitoring  of  a 
significant  number  of  target  events. 

2.  TECHNICAL  APPROACH 

FARM  Database 

Our  initial  efforts  concentrated  on  assembling  a  database  of  wavefonns  from  the  IRIS 
FARM  archive,  which  consists  of  broadband  data  from  the  global  seismic  networks  as 
well  as  portable  seismic  arrays  deployed  in  PASSCAL  experiments.  This  involves 
transferring  the  data  from  the  IRIS  DMC  and  converting  waveforms  to  SAC.  We  then  run 
programs  for  quality  assurance,  redundancy  checks,  instrument  response  corrections, 
meta-data  archiving,  and  signal  rotation  to  tangential  and  radial.  The  SAC  files  are 
currently  stored  on  a  RAID  system,  with  backup  for  redundancy  in  the  event  of  hardware 
failures. 

We  have  also  computed  index  files  and  signal-to-noise  estimates  for  the  P  and  S 
arrivals,  which  facilitates  later  processing  and  also  provides  a  check  on  the  timing 
integrity  of  the  SAC  waveforms.  Our  basic  approach  is  to  measure  the  signal-to-noise  as 
the  ratio  between  the  maximum  amplitude  (peak  to  trough)  in  a  time  window  that 
contains  the  phase  of  interest  and  a  pre-event  noise  window  of  equal  length.  Because  the 
raw  data  are  broadband,  we  perfonn  this  operation  separately  for  different  frequency 
bands. 
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Event  Term  v.  Distance 


Figure  3.  Distance-dependent  Earth-response  terms  for  P- wave  amplitude  from  10  to  90 
degrees. 

Amplitude  Analysis 

Triplicated  phases  are  difficult  to  distinguish  from  each  other,  overlapping  in  time, 
and  having  varying  amplitudes.  This  is  further  compounded  by  unconstrained  source¬ 
time  functions  of  the  earthquake,  and  heterogeneous  structure  beneath  seismic  stations 
causing  waveform  variability.  Here,  we  attempt  to  remove  the  source  and  station  tenns 
from  all  measured  amplitudes  for  high  signal-to-noise  P  waves,  to  constrain  the 
amplitude  variation  that  results  from  seismic  structure.  Initially  we  assume  as  in  Figure  4 
that  the  maximum  peak-to-peak  log  amplitude  (max-min),  a„  of  the  Ith  P- wave  is  the  sum 
of  a  source  tenn,  sh  a  receiver  term,  rk,  and  an  Earth  response  tenn,  ei. 


a,  =  Sj  +  rk  +  <?/. 


Given  Nu  E-wave  amplitudes,  from  any  Ns  sources,  recorded  at  Nr  sources,  it  is 
possible  to  solve  for  Ne  Earth-response  terms  as  a  function  of  distance.  The  initial  source 
tenn  for  each  event  is  the  mean  of  all  at  values  for  that  event.  The  initial  receiver  terms 
are  the  means  of  all  amplitude  measurements  for  that  receiver  minus  the  corresponding 
source  terms  (rk  =  a,  -  Sj ).  The  distance-dependent  Earth-response  tenn  is  the  mean  of  all 
amplitude  measurements  for  each  distance  minus  the  corresponding  source  and  receiver 
tenns  (ei  =  a,  -  sj  -  n ).  We  can  then  iteratively  solve  for  more  accurate  source,  receiver, 
and  Earth  response  tenns,  with 

sj  =  at  -  rk  -  ei 
rk  =  at  -  sj  -  ei 
ei  =  at  -  sj  -  rk. 
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Figure  4.  A  graphical  representation  of  the  sources  of  seismic  signals  (event,  Earth 
response  term,  receiver  term,  and  noise)  that  we  attempt  to  decompose. 
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Figure  5.  Examples  of  three  different  types  of  perturbations  to  transition  zone  structure: 
(a)  transition  zone  thickness,  (b)  velocity  contrast,  and  (c)  interface  thickness.  The 
corresponding  WKB J  maximum  amplitudes  for  P  waves  (d-f)  change  as  a  function  of 
distance.  The  models  are  shown  as  black,  blue,  and  red  lines.  The  data  are  the  dashed 
lines. 
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When  applied  to  all  P-wave  data  in  the  FARM  database,  this  method  provides  an 
estimate  of  amplitude  as  a  function  of  distance,  as  shown  in  Figure  3.  The  overall 
amplitude  versus  distance  behavior  is  very  similar  to  that  obtained  by  Veith  and  Clawson 
(1972)  for  magnitude  calibration.  The  13°  to  33°  distance  range  stands  out  as  anomalous 
compared  to  the  trend  of  decreasing  amplitude  with  distance. 

We  have  adapted  a  WKBJ  synthetic  seismogram  code  (Chapman,  1978)  to  generate 
suites  of  model  predictions  for  comparisons  to  the  data.  Our  analyses  have  concentrated 
on  characterizing  the  globally  averaged  wavefield  as  a  starting  reference  point  for  studies 
of  regional  variations.  By  comparing  the  Earth-response  amplitude  terms  to  WKBJ 
synthetic  waveform  amplitudes,  it  is  possible  to  detennine  which  seismic  models  are 
most  probable.  The  amplitudes  in  the  13°  to  33°  range  are  sensitive  to  the  seismic 
velocity  structure  of  the  transition  zone.  Amplitudes  are  highly  sensitive  to  interface 
depth,  sharpness,  and  velocity  contrast  for  each  of  the  discontinuities  {410,  520,  and 
660).  Figure  5  shows  how  amplitude  can  vary  with  small  changes  in  seismic  velocity 
structure.  While  not  modeled  here,  the  amplitudes  are  also  strongly  dependent  upon 
attenuation  structure. 


3.  RESULTS  AND  DISCUSSION 

We  model  the  data  with  a  mass-forward  modeling  algorithm  called  the  Niching 
Genetic  Algorithm  (NGA)  {Koper  et  al.,  1999).  The  NGA  uses  an  evolutionary  paradigm 
to  create,  mutate,  and  breed  a  population  of  velocity  structure  models  such  that  those  with 
poor  fits  to  the  data  are  eliminated,  and  those  with  good  fits  survive  to  continue  the 
search  of  the  model  space  for  the  best  solution.  Using  this  paradigm,  we  find  the  optimal 
solution  to  the  P-wave  data  (Figure  6).  Using  the  NGA,  we  parameterize  velocity 
contrast,  density  contrast,  interface  depths,  and  interface  thickness  for  each  interface 
(410,  520,  and  660)  as  shown  in  Figure  6.  The  sub-660  gradient  is  defined  by  the 
interface  at  the  bottom  of  the  660  and  a  peg  depth  parameter  below  which  the  P-wave 
velocity  and  density  are  equal  to  AK135. 

The  resultant  model  (Figure  7)  is  similar  to  the  P- wave  velocity  structure  found  by 
Lawrence  and  Shearer  (2006)  using  NGA  waveform  modeling  of  reflected  energy  from 
the  transition  zone  (PdP,  SdS,  Pds,  and  Ppdp).  Both  models  have  roughly  the  same 
velocity  contrasts  at  the  410,  and  660  km  discontinuities.  The  sub-660  gradient  is 
identical  within  the  resolution  of  each  study.  The  estimated  amplitude  and  depth  of  the 
520  km  discontinuity  are  more  different,  but  are  poorly  constrained  in  both  studies 
relative  to  the  larger  interfaces.  In  fact,  the  synthetic  Earth-response  term  of  the  model 
from  Lawrence  and  Shearer  (2006)  provides  an  equivalent  fit  to  the  observed  global 
Earth  response  as  the  most  optimal  model  of  this  study. 
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Figure  6.  Parameterization  of  the  P  velocity  model.  Each  interface  (410,  520,  &  660)  can 
vary  in  depth,  thickness,  and  velocity  contrast. 


VP  (km/s) 

Figure  7.  The  (red)  best  fit  model  after  100  generations  of  an  NGA  simulated  inversion 
yields  a  similar  result  to  a  similar  parameter  search  based  on  (blue)  PdP,  SdS,  Pds,  and 
Ppdp  reflected  phases  (Lawrence  and  Shearer,  2006). 
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Wavefield  Analysis 

In  the  amplitude  section,  we  neglected  time  other  than  to  ensure  that  we  have 
windowed  the  appropriate  phase.  Here,  we  introduce  a  similar  analysis  that  iteratively 
solves  for  time-dependent  source,  station,  and  Earth-response  terms.  Mathematically,  one 
can  conceptualize  a  waveform,  a{t),  as  the  convolution  of  the  source  term,  s(t),  the 
receiver  term,  r(t),  and  the  distance-dependent  Earth-response  term,  e(D ,t),  plus  noise, 
n(t). 


a(?)  =  s(?) *  r(7) *  e(A, ?)+  «(?) 

Here,  we  assume  that  the  source  term,  s,  is  a  combination  of  the  earthquake  rupture  and 
the  near-source  Earth  response  function.  Furthermore,  the  receiver  term  is  idealized  as  the 
convolution  of  the  instrument  response  and  the  near-receiver  Earth  response. 

Ideally,  one  could  estimate  the  deep-Earth  response  by  stacking  waveforms  with  the 
source  and  receiver  terms  (s(t)  and  r(t))  deconvolved.  Unfortunately,  the  source  and 
receiver  terms  are  poorly  constrained,  the  waveforms  are  not  easily  stackable  because  of 
dissimilarity  and  polarity  issues,  and  noise  is  typically  problematic.  Therefore  we  develop 
and  employ  an  automated  stepwise  stacking  and  deconvolution  technique  to  boost  the 
signal-to-noise  ratio  and  to  solve  for  the  time-dependent  source  and  receiver  terms.  We 
improve  waveform  similarity  by  examining  the  instantaneous  amplitude  or  envelope 
function  of  the  wavefield,  ai{t),  rather  than  the  raw  waveform.  We  initially  estimate  the 
source,  receiver,  and  Earth-response  functions  as  delta  functions  with  amplitudes 
provided  by  the  amplitude  analysis  above. 

With  preliminary  estimates  of  s,  r,  and  e,  we  can  iteratively  solve  for  each  term  from 
the  envelope  function, 


^(0  =  r(t)*'1  e(A,t) 

r{t)  =  aE(')*~'  s(f)*'1  e(A,t) 
e(A,t)  =  aE(t)*-'  r(t ) 

With  this  analysis,  we  obtain  a  time-  and  distance-dependent  Earth-response  term. 

For  all  data  in  the  FARM  database,  the  Earth-response  term  is  shown  in  Figure  8  as 
determined  with  spectral  deconvolution.  In  addition  to  spectral  division  techniques,  we 
investigate  the  benefits  of  employing  iterative  deconvolution.  While  spectral 
deconvolution  is  computationally  faster  than  iterative  deconvolution,  it  is  acausal,  often 
noisy,  and  results  depend  heavily  on  a  waterlevel.  Iterative  deconvolution  is  slower,  but  is 
causal,  can  solve  for  positive  peaks  only,  and  provides  a  minimum  information  solution. 
In  Figures  8  and  9  we  compare  results  of  the  spectral  division  and  iterative  deconvolution 
techniques. 
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Figure  8.  The  time-  and  distance-dependent  Earth-response  function  calculated  from  all 
FARM  F’-wavc  data  using  the  spectral  deconvolution  method.  The  box  outlines  the 
region  of  triplicated  data.  Blue  is  positive  amplitude,  red  is  negative.  Various  phases  are 
indicated  as  P,  PP,  and  PcP. 

There  appear  to  be  significant  differences  between  the  two,  but  the  majority  of  the 
differences  are  due  to  the  longer-period  smoothing  necessary  for  the  spectral  division 
technique.  Consequently  some  features  appear  to  be  sharper  in  the  iterative  deconvolution 
technique.  For  example,  the  iterative  deconvolution  technique  provides  a  sharper  high 
amplitude  anomaly  at  25  to  28  degrees  after  the  initial  pulse  (3  to  5  seconds).  This  pulse 
is  associated  with  the  P  branch  of  the  660  km  discontinuity. 
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Figure  9.  The  time-  and  distance-dependent  Earth-response  function  calculated  from  all 
FARM  T’-wavc  data  using  the  iterative  deconvolution  method.  The  box  outlines  the 
region  of  triplicated  data.  Blue  is  positive  amplitude,  red  is  negative.  Various  phases  are 
indicated  as  P,  PP,  and  PcP. 

While  the  iterative  deconvolution  technique  provides  resolution  at  higher  frequencies, 
the  differences  between  the  modeled  synthetics  for  each  technique  shows  that  there  is 
minimal  improvement  in  resolution  relative  to  the  spectral  method.  While  all  further 
analyses  were  conducted  with  both  the  spectral  division  technique  and  the  iterative 
deconvolution  technique,  we  only  present  results  from  one  method  because  of  their  great 
overall  similarity. 

The  amplitudes  of  the  globally  stacked  Earth  response  function  using  iterative 
deconvolution  at  time  zero  for  each  distance  grossly  match  the  estimates  from  the  simpler 
maximum  amplitude  decomposition  (see  Figure  10).  Both  the  raw  amplitudes  and  the 
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Earth-response  term  amplitudes  have  elevated  amplitudes  between  10  and  30  degrees.  Of 
particular  interest,  the  low  amplitudes  observed  at  ~30  degrees  in  both  amplitude 
distributions  is  continued  to  nearly  40  degrees  in  zero-time  Earth-response  term 
amplitudes.  At  ~40  degrees,  there  is  a  small  elevated  amplitude.  Variations  at  ~35 
degrees  result  from  amplitude  differences  in  the  first  arrival  relative  to  the  later  arriving  P 
branch  from  the  660  km  discontinuity.  The  amplitudes  gradually  decrease  from  ~40 
degrees  to  90  degrees.  Beyond  90  degrees  the  amplitudes  decrease  rapidly  because  of  the 
influence  of  the  low  velocities  in  the  outer  core. 
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Figure  10:  The  amplitudes  of  the  globally  stacked  Earth  response  function  using  iterative 
deconvolution  at  time  zero  for  each  distance  (yellow)  grossly  match  the  average  estimates 
from  simple  maximum  amplitudes  for  P- waves  (red).  The  black  dots  illustrate  the  scatter 
in  the  P- wave  amplitudes. 
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Figure  1 1 .  This  figure  shows  (a)  the  synthetic  Earth-response  function  for  the  seismic 
velocity  model  found  by  Lawrence  and  Shearer  (2006)  (LS06),  and  deviations  in  Earth- 
response  functions  relative  to  LS06  as  a  result  of  (b)  a  deeper  660  interface,  (c)  a  larger 
Fp660  contrast,  and  (d)  a  thick  410  interface.  Blue  is  positive,  red  is  negative. 

To  model  these  results,  we  calculate  synthetic  time-  and  distance-dependent  Earth- 
response  functions  using  WKBJ.  To  estimate  the  appropriate  frequency  content  and 
heterogeneity,  we  deconvolve  the  E-wave  slant-stacked  synthetic  waveforms  from  45  to 
80  degrees  and  convolve  the  stacked  observed  waveforms  for  the  same  distance  range. 
Synthetics  using  the  model  of  Lawrence  and  Shearer,  (2006)  (LS06)  are  similar,  but  not 
identical  to  the  observed  Earth-response  functions  (Figure  1  l.a).  Variants  from  LS06 
demonstrate  how  the  Earth-response  function  changes  as  a  result  of  structure  (Figure 
ll.b-d). 
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Figure  12:  A  comparison  showing  that  the  Earth-response  function  calculated  from  the 
model  of  Lawrence  and  Shearer,  (2006)  is  very  similar  to  the  global  Earth  response 
function.  These  figures  were  calculated  with  decomposition  using  spectral  division. 

As  shown  in  Figure  12,  the  optimal  P-wave  model  found  with  a  niching  genetic 
algorithm  simulated  inversion  for  the  wave  Held  is  remarkably  similar  to  the  most  optimal 
model  found  in  Figure  7.  Note  that  the  most  optimal  model  is  merely  representative  as  we 
assume  that  3D  heterogeneity  has  a  broadening  effect,  and  that  we  do  not  know  precisely 
the  scale  of  this  heterogeneity.  The  primary  difference  between  the  optimal  model  found 
with  wavefield  analysis  and  that  from  amplitude  analysis  is  that  the  estimated  depths  for 
the  global  model  are  slightly  different,  which  accounts  for  the  timing.  Recall  that  timing 
is  inherently  included  in  the  wavefield  analysis  and  not  in  the  amplitude  analysis.  We 
expect  that  this  method  will  also  work  with  shear  waves,  although  so  far  we  have  only 
examined  P  waves. 

Event  Characterization 

Once  the  station  terms  are  well  defined  with  the  global  analysis,  it  is  possible  to 
perform  regional  amplitude  or  wavefield  analysis  for  high-quality  individual  events  or 
sets  of  events.  With  less  data  it  is  critical  to  use  only  high-quality  events  and  well- 
calibrated  stations. 

By  deconvolving  the  Earth  response  term  and  station  response  term  from  the 
envelopes  of  data  from  a  single  event,  it  is  possible  to  determine  the  magnitude  based  on 
scaling  to  the  observed  event  terms  (Figure  13).  By  scaling  the  global  Earth  response 
term’s  zero-time  amplitudes  to  the  observed  zero-time  amplitudes,  it  may  be  possible  to 
estimate  the  event  magnitude  using  distances  between  the  13  to  33  degree  distance  range. 
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Figure  13:  The  Earth  response  term  for  a  single  event,  as  (left)  traces,  and  (right)  as  the 
zero-time  amplitude  (x’s)  as  a  function  of  event- to-station  distance 

The  scaling  factors  provide  an  estimate  of  earthquake  magnitude.  The  scaling  factor 
(sf)  is  correlated  with  body  wave  magnitude  (Mb)  with  an  R  value  of  0.7  for  9  of  the  12 
tested  events  in  2005  as  shown  in  Figure  14.  Because  the  Earth-response  term  is 
calculated  with  events  from  1990  to  2004,  the  agreement  between  Mb  and  sf  is  not  a 
result  of  pre-constrained  relations  for  these  events.  Mb  and  sf  are  related  with  the 
equation  of  sf=  1 .28Mb  -  8.2.  Note  that  the  results  of  3  events  failed  due  to  lack  of 
agreement  (correlation)  between  the  event  and  global  Earth-response  functions.  These 
failures  stemmed  from  more  complex  sources.  The  9  more  impulsive  sources  provided 
reasonable  correlations  between  the  new  Earth-response  tenn  and  the  global  Earth 
response  term. 

By  using  regional  Earth-response  functions  the  event  magnitude  estimation  may 
become  more  accurate.  The  regional  Earth-response  tenns  will  likely  be  more 
representative  of  regional  structure  and  include  less  damping/broadening  effects  due  to 
heterogeneity. 
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Figure  14:  The  relative  scaling  factor  between  the  zero-time  Earth-response  amplitudes  of 
a  single  event  and  the  zero-time  Earth-response  amplitudes  of  the  global  dataset  plotted 
versus  body  wave  magnitude. 

4.  CONCLUSIONS 

This  analysis  has  shown  that  it  is  possible  to  divide  amplitude  terms  from  digital 
seismic  records  into  source,  receiver,  and  Earth-response  terms  that  account  for  the 
complicated  wavefield  variations  in  the  13°  and  33°  distance  range.  Through  a  multi¬ 
dimensional  component  analysis,  each  term  is  obtained  iteratively.  With  well-resolved 
receiver  terms  and  Earth-structure,  it  may  be  possible  to  estimate  source  functions  with 
data  at  these  regional  distances.  Note  that  the  method  requires  large  numbers  of 
seismograms  from  many  stations  to  achieve  reliable  results.  It  is  only  in  the  past  few 
years  that  global  data  and  station  coverage  have  increased  to  a  sufficient  level  that  such 
an  analysis  can  be  roughly  conducted.  In  the  future,  improved  results  are  likely  as  more 
data  are  collected  from  regions  of  increased  station  coverage. 

By  using  a  priori  constraints  on  the  Earth-response  term  and  station  terms  from  prior 
event  databases,  it  is  possible  to  detennine  the  rough  magnitude  of  some  new  events 
using  the  13-33  degree  distance  amplitudes.  This  method  does  not  work  well  for  all  new 
events  because  of  complex  sources,  but  works  well  for  more  impulsive  sources. 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


DMC 

Data  Management  Center 

FARM 

Fast  Archive  Recovery  Method 

IRIS 

Incorporated  Research  Institutions  in  Seismology 

NGA 

Niching  Genetic  Algorithm 

PASSCAL 

Program  for  the  Array  Studies  of  the  Continental  Lithosphere 

RAID 

Redundant  Array  of  Inexpensive  Disks 

SAC 

Seismic  Analysis  Code 

WKBJ 

W  entzel-Kramers-Brillouin- J  effreys 
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